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In this work we put forward an exact one-particle framework to study nano-scale Josephson junc- 
tions out of equilibrium and propose a propagation scheme to calculate the time-dependent current in 
response to an external applied bias. Using a discrete basis set and Peierls phases for the electromag- 
netic field we prove that the current and pairing densities in a superconducting system of interacting 
electrons can be reproduced in a non-interacting Kohn-Sham (KS) system under the influence of 
different Peierls phases and of a pairing field. In the special case of normal systems our result pro- 
vides a formulation of time-dependent current density functional theory in tight-binding models. An 
extended Keldysh formalism for the non-equilibrium Nambu-Green's function (NEGF) is then intro- 
duced to calculate the short- and long-time response of the KS system. The equivalence between the 
NEGF approach and a combination of the static and time-dependent Bogoliubov-deGennes (BdG) 
equations is shown. For systems consisting of a finite region coupled to Af superconducting semi- 
infinite leads we numerically solve the static BdG equations with a generalized wave-guide approach 
and their time-dependent version with an embedded Crank-Nicholson scheme. To demonstrate the 
feasibility of the propagation scheme we study two paradigmatic models, the single-level quantum 
dot and a tight-binding chain, under dc, ac and pulse biases. We provide a time-dependent picture 
of single and multiple Andreev reflections, show that Andreev bound states can be exploited to 
generate a zero-bias ac current of tunable frequency, and find a long-living resonant effect induced 
by microwave irradiation of appropriate frequency. 
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I. INTRODUCTION 



In the last two decades superconducting nanoelectron- 
ics has emerged as an interdisciplinary field bridging dif- 
ferent areas of physics like superconductivity, quantum 
transport and quantum computation^— For practical ap- 
plications the reduction of heat losses in superconducting 
circuits constitutes a major advantage over semiconduc- 
tor electronics where a molecular junction is more subject 
to thermal instabilities^— 

The idea of exploiting atomic-size quantum point con- 
tacts or quantum dots coupled to superconducting leads 
as quantum bits (QUBIT) has received significant atten- 
tion both theoretically and experimentally^— The state 
of a QUBIT evolves in time according to the Schrodingcr 
equation for open quantum systems and can be manip- 
ulated using electromagnetic pulses of the duration of 
few nano-seconds or even faster. Due to the reduced 
dimensionality and the high speed of the pulses these 
systems can be classified as ultrafast Josephson nano- 
junctions (UF-JNJ). The microscopic description of the 
out-of-equilibrium properties of an UF-JNJ is not only 
of importance for their potential applications in future 
electronics but also of considerable fundamental interest. 
The quantum nature of the nanoscale device leads to a 
sub-harmonic gap structure jiSr— ac characteristics; 1 ^ 18 
current-phase relation ; 19 i 20 etc. that differ substantially 
from those of a macroscopic Josephson junction. Furt her- 



more, there are regimes in which the electron-electron 
scattering inside the device plays an important rolejSir— 
We here focus on a different relevant aspect of UF- 
JNJ, namely the ab initio description of their short time 
responses. Considerable theoretical progresses have been 
made to construct a first-principle scheme of electron 
transport through molecules placed between normal met- 
als. On the contrary, despite the recent experimental 
advances in fabricating superconducting quantum point 
contacts, a first-principle approach to superconducting 
nanoelectronics is still missing. Furthermore, time- 
dependent (TD) properties like the switch on/off time of 
the current or the response to time-dependent ac fields 
or train pulses has remained largely unexplored. There 
are several difficulties related to the construction of a 
feasible time-dependent approach already at a mean-field 
level. The system is open, the electronic energy scales are 
2-3 orders of magnitude larger than a typical supercon- 
ducting gap, the problem is intrinsically time-dependent 
(even for dc biases), and the possible formation of An- 
dreev bound states (ABS) give rise to persistent oscilla- 
tions in the density and current. The time-evolution of lo- 
calized wave-packets scattering across a superconductor- 
normal interface was explored long agOi2£r— More re- 
cently the analysis has been extended to scattering states 
in superconductor-device-normal (S-D-N) junctions us- 
ing the wide-band-limit (WBL) approximation^ and in 
superconductor-device-superconductor (S-D-S) junctions 
by approximating the leads with finite size reservoirs^ 
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However, there has been no attempt to calculate the re- 
sponse of S-D-S junctions to TD applied voltages using 
truly semi-infinite leads. 

In this work we propose a one-particle framework to 
study TD quantum transport in UF-JNJ, construct a 
suitable propagation scheme and apply it to study gen- 
uine TD properties like the switch on/off of the current, 
the onset of a Josephson regime, ABS oscillations, ac 
transport and the time-evolution of multiple Andreev re- 
flections. 

The one-particle framework, described in Section III Al 
and III Bt is an extension of TD superconducting density 
functional theory^ to systems with a discrete basis and 
is built on the mapping from densities to potentials pro- 
posed by van Leeuwen^ and Vignale^ It is shown that 
under reasonable assumptions the current density and 
pairing density of an interacting system perturbed by a 
TD electromagnetic field can be reproduced in a Kohn- 
Sham system of non-interacting electrons perturbed by 
a TD electromagnetic and pairing fields, and that these 
fields are unique. In the special case of normal systems 
such result provides a formulation of TD current density 
functional theory in tight-binding models. 

An extended Kcldysh formalism for the non- 
equilibrium Nambu-Green's function is introduced in 
Section III CI and used to calculate the time-dependent 
current, density and pairing density of the Kohn-Sham 
Hamiltonian. By adding a vertical imaginary track to 
the original Keldysh contour— ~— we are able to extract 
the response of the system just after the application of 
the bias (transient regime) and to describe the onset of 
the Josephson regime. We also show the equivalence be- 
tween the equations of motion for the Nambu-Grccn's 
function on the extended contour and the combination 
of the static and TD Bogoliubov-DcGennes equations. 

In Section ITTTI we illustrate a procedure for the calcula- 
tion of the one-particle eigenstates of a system consisting 
of M semi-infinite superconducting leads coupled to a fi- 
nite region C. These states are then propagated in time 
according to the TD Bogoliubov-DcGennes equations us- 
ing an embedded Crank-Nicholson algorithm which re- 
duces to that of Refs. I37ll38l in the case of normal leads. 
The propagation scheme is unitary (norm conserving) 
and incorporates exactly the transparent boundary con- 
ditions. 

The feasibility of the method is demonstrated in Sec- 
tion [IV] where we calculate the TD current, density and 
pairing density of S-D-S junctions under dc, ac and pulse 
biases. The paradigmatic model with a single atomic 
level connected to a left and right superconducting leads 
is investigated in detail. We provide a time-dependent 
picture of single and multiple Andreev reflections and of 
the consequent formation of Cooper pairs at the inter- 
face. We show that the smaller is the bias the longer and 
the more complex is the transient regime. We also study 
how the system relaxes after the bias is switched off. Due 
to the presence of ABS a tiny difference in the switch-off 
time can cause a large difference in the relaxation be- 



havior with persistent oscillations of tunable frequency. 
ABS also play a crucial role in microwave ac transport. 
Tuning the frequency of the microwave field according 
to the ABS energy difference one produces a long-living 
transient resonant effect in which the amplitude of the 
ac current is about an order of magnitude larger than 
that of the current out of resonance. Finally we consider 
one-dimensional atomic chains coupled to superconduct- 
ing leads. We calculate the TD current density pattern 
along the chain for dc (ac) biases and show a clear-cut 
transient scenario of the multiple (photon-assisted) An- 
dreev reflections. A summary of the main findings and 
an outlook on future perspectives are drawn in Section 

E 

II. GENERAL FORMULATION 

A. Hamiltonian of the system 

The Hamiltonian of a system of interacting electrons 
can be written in terms of the field operators Vv(r) 
(V>J.(r)) which destroy (create) an electron of spin a 
in position r. We expand the field operators in some 
suitable basis of localized orbitals y m (r) as Vv(r) = 
Em'iMVmir)' Assuming, for simplicity, that the (fim's 
are orthonormal the c's operators obey the anticommu- 
tation relations 

{^mcr;C na ./} — o'aa'O'nm- (1) 

In the presence of an external static electromagnetic and 
pairing field the Hamiltonian has the general form 

i?o = #o + Ao + AJ + F int . (2) 

The first term is the free-electron part and reads 

Ko = T '»« el7m " d iiaCna (3) 

a ran 

with real symmetric hopping parameters T mn = T nm and 
real antisymmetric phases j mn = — j nm . The phases ac- 
count for the presence of an external vector potential 
A(r), in accordance with the Peierls prescription. If we 
use a grid basis for the expansion of the field operators 
with grid points r m then "f mn — ^ J r m dl- A(r). The sec- 
ond term in Eq. @ represents the pairing field operator 
which couples the pairing density operator to an external 
field and reads 

A = ^A m c] rit c I t ni . (4) 

III: 

We notice that the pairing field A m is local in the chosen 
basis. This term is usually set to zero since the transition 
to a superconducting state is caused by the interaction 
part. Our motivation to include it at this stage will soon 
become clear. The interaction part of the Hamiltonian 
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Hint contains terms more than quadratic in the c's op- 
erators. We do not specify the form of H lnt which can 
be any. We, however, require that it commutes with the 
density operator h ma = & ma c mG 

[.Hint, h ma ] = 0, Vm,cr. (5) 

The above condition is fulfilled on a grid basis as well as 
in tight-binding models with Hubbard-like interactions. 

We are interested in the dynamics of the system when 
an extra time-dependent electromagnetic field and pair- 
ing potential is switched on at t — 0. The pairing po- 
tential must here be considered as an independent ex- 
ternal field. Since the time-dependent part of the scalar 
potential can always be gauged away we restrict to time- 
dependent Hamiltonians of the form 

H(t) = K(t) + A(t) + A*(t) + H in t, (6) 

where 

k ® = E E T mn e l ^cl a c n<J (7) 

ex run 

and 

A(t)=X)A w (t)4,t44- ( g ) 

tn 

In 1994 Wackcr, Kiimmel and Gross^I put forward a 
rigorous framework, known as TD Density Functional 
Theory for Superconductors (SCDFT), to study the dy- 
namics of a superconducting system in the continuum 
case. The continuum Hamiltonian can be obtained from 
the Hamiltonian in Eq. ((B) with the tp m 's a grid ba- 
sis in the limit of zero spacing. They proved that given 
an initial many-body state |$o) the current and pair- 
ing densities evolving under the influence of two different 
vector potentials A and A' and/or two different pair- 
ing potentials A and A' are always different. This re- 
sult renders all observable quantities functionals of the 
current and pairing densities, which can therefore be cal- 
culated in a one-particle manner^ The original formu- 
lation relies on the assumption that the time-dependent 
current and pairing densities of the interacting Hamil- 
tonian can be reproduced in a non-interacting Hamil- 
tonian under the influence of another vector and pair- 
ing potential, i.e., that the interacting A- A densities are 
also non-interacting A-A representable. The interact- 
ing versus non-interacting representability assumption is 
present also in the original formulation of TD Density 
Functional Theory (DFT) by Runge and Gross^ and TD 
Current Density Functional Theory (CDFT) by Ghosh 
and Dhara4£ The representability problem in TDDFT 
was solved by van Leeuwen who proved that the TD den- 
sity of a system with interaction H lnt under the influence 
of a TD scalar potential V can be reproduced in another 
system with interaction H[ nt under the influence of a TD 
scalar potential V' and that V' is unique.— We will re- 
fer to such result as the van Leeuwen theorem. Taking 



H[ nt — the van Leeuwen theorem implies that the TD 
interacting density can be reproduced in a system of non- 
interacting electrons. Later Vignale extended the van 
Leeuwen theorem to solve the representability problem in 
TDCDFT<22 In the next section we show that the results 
by van Leeuwen and Vignale can be further extended 
to solve the representability problem in TDSCDFT. The 
theory is formulated on a discete basis and it is not lim- 
ited to pure states, implying that we also have access to 
the finite-temperature domain. 

B. The one-particle Kohn-Sham scheme of 
TDSCDFT 

Let p(f) be the density matrix at time t of the system 
described by the Hamiltonian in Eq. (JB|). We denote 
by 0(t) = Tr {p(t)0(t)} the time-dependent ensemble 
average of a generic operator 0(t), where the "Tr" sym- 
bol signifies the trace over a complete set of many-body 
states. The average 0(t) obeys the equation of motion 

jO{t) = ^0(t) + iTr{p(t)m), 0(t)}}. (9) 

It is easy to verify that when 0(t) is the density operator 
n™, = J2a 4w4<ri Eq. © yields 

j t n m (t) = J2 J "»> W - 4Im [A; n (i)P m We- 2lT — '] , 

n 

(10) 

where J mn (t) and P m (f) are the expectation value of the 
bond-current operator 

Jmn(f) = ~ ^ ^ {Pmn^-^ ^ ^ C^^CrKj H.C^ (I I ) 

cr 

and pairing density operator 

P m (t) = c rni c mt e 2t £ dt ' T — = c ml c mf e 2iT ™™ t . (12) 

Equation (|I0j) is the proper extension of the continuity 
equation to systems exposed to a pairing field. The term 
A(i) + AT(t) acts as if there were TD sources and sinks. 

Notice that under the gauge transformation c na — > 
e i/3„(t)g nff ( w ith /3„(0) = 0) the on-site energies change as 
T mm — > T mm — df3 m (t)/dt while the phases and the pair- 
ing field change according to 7„ m (i) -> Jmn(t) + Pm(t) - 
P n (t) and A m (t) -> A m (t) exp[2i0 m (t)]. Therefore the 
bond-current operator J mn and pairing density operator 
P m are gauge invariant. In a grid basis representation 
with grid points r m the phases /3 rn (t) are the discretized 
values of the scalar function A(r m ,t) which defines the 
gauge-transformed vector potential A and scalar poten- 
tial V: A -> A + cVA and V -> V - dA/dt. 

The equation of motion for the bond-current J mn (i) 
can be cast as follows 

■j;J mn (t) = K mn (t)-^-j mn (t) + F mn (t). (13) 



4 



The first term in the r.h.s. is exactly dJ mn (t)/dt; the 
operator K mn (t) = J2a ( r m n e l7m " (t) cJ n(T c„ CT + H.c.) is 
the energy density of the bond m-n. The second term 
in the r.h.s. is, therefore, the average of F mn (t) = 
i[H{t),J mn (t)\, sec Eq. ©. 

The derivation of the equation of motion for the pairing 
density P m (t) is also straightforward and leads to 

(j t - 2iT mm j P m (t) = iA m (t)[n m (t) - l] e 2lT — * 

+ iG m (t)e 2lT ™ *, (14) 

with G m (t) = [K(t) + H int , Cml&mf]. 

We now ask the question whether the densities J mn (t) 
for all bonds m-n with T mn ^ and P m {t) can be re- 
produced in a system with a different interaction Hamil- 
tonian H' int under the influence of TD phases ^'{t) and 
pairing potential A'(t) starting from an initial density 
matrix /S'(0). 

For the densities to be the same at time t = we have 
to choose p'(0) and 7'(0) in such a way that 

Tr{p'(0)j; m (0)} = Tr{p(0)J m „(0)}, (15) 

Tr{p'(0)P m (0)} = Tr{p(0)P m (0)}. (16) 

Notice that in the primed system the bond-current op- 
erator J' mn is different from J mn since the phases 7' are 
generally different from 7. On the contrary the pairing 
density operator is the same in the two systems. Equa- 
tions (|15ll6p define the compatible initial configurations 
of the primed system. 

We answer the above question affirmatively by showing 
that given a compatible initial configuration [p'(0), 7'(0)] 
and under reasonable conditions there exist "f'(t) and 
A'(t) for which the bond-current and pairing density of 
the original and primed system are the same at all times. 
The formal statement is enunciated in the following 

Theorem : Given a compatible initial configuration 
[p'(0),7'(0)] such that 

K' mn (0) = T?{p'{Q)Y^{T mn e l ^)cl a c na +H.c.)} ^ 

(17) 

for all bonds m-n with T mn ^ 0, and 

n' m (0)=Tr{p'(0)h m }^l, (18) 

which implies that at time t = none of the orbitals ip m 
are half filled in the primed system, there exist a unique 
set of continuous phases 7'(t) and pairing potential A'(t) 
that reproduce in the primed system the densities J mn (t) 
and P m (t) of the original system. 

Remarks : Before presenting the proof of the The- 
orem we discuss few relevant implications. (1) If the 
original system is a superconducting system with an at- 
tractive interaction H lnt and a vanishing pairing field, 
i.e., A = 0, the theorem implies that the bond-currents 



and pairing densities can be reproduced in a system of 
non- interacting electrons, i.e., H' lnt = perturbed by 
TD phases 7' and pairing field A'. In the following we 
will refer to such non-interacting system as the Kohn- 
Sham (KS) system and to the TD perturbation as the 
KS phases and KS pairing potential. In Section IIIII 
we describe how to perform the time-evolution of such 
KS systems for geometries relevant to quantum trans- 
port. (2) For interacting systems with A = and ini- 
tially in equilibrium in the absence of electromagnetic 
fields the phases 7(0) = and hence J mn (0) — for all 
bonds. In the KS system a possible compatible initial 
configuration is therefore 7'(0) = and p'(0) such that 
the expectation value of the one-particle density matrix 
n mn(0) = S CT ^{p'WcLa^na} is real. For such initial 
configurations the condition (|T7|) becomes n' mn (0) 7^ 
for all bonds m-n with T mn 7^ 0. (3) If we ask the ques- 
tion whether only the bond-currents J m n(t) of a system 
with Hamiltonian ((6j> and zero pairing field, i.e., A = 0, 
can be reproduced in a system with zero pairing field, 
i.e., A' = 0, and different interactions H[ nt under the in- 
fluence of different phases 7' starting from some initial 
density matrix p'(0), the answer is affirmative provided 
that p'(0) and Y(0) fulfill Eqs. (|15I17|) . This corollary ex- 
tends TDCDFT to tight-binding models using the Peierls 
phases as the basic KS fields and lays down the basis for 
a density functional TD theory in discrete systems^ 

We conclude this Section with the proof of the Theo- 
rem. 

Proof : The current and pairing densities of the 
primed system obey the equations of motion (|13I14[) with 
K mn (t) K' mn (t), F mn (t) Fl nn (t) and n m (t) 
n' m (t), G m (t) — > G' m {t). Therefore, for a generic time 
t the densities of the two systems are the same provided 
that 

^mnW^K*) = K mn (t)—-f mn (t) 

+ F mn {t)-F: nn {t), (19) 



[n' m (t) - l]A' m {t) = [n m (t) - l]A m (i) 

+ G m {t)-G' m {t). (20) 

A discussion on the existence and uniqueness of the so- 
lution for the coupled Eqs. (|19H20[) is rather complicated 
since the dependence on the phases 7' and potentials A' 
in F' and G' enters implicitly via the TD density matrix 
p'(t). To proceed further we then follow the approach 
of Vignale and assume that the time-dependent phases 
and pairing potentials and hence all expectation values 
are analytic functions of time around t = Expand- 
ing all quantities in Eqs. (|19H20I) in their Taylor series 
and equating the coefficients with the same power of t 
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we obtain 

l-i 

(Z + l)ifmn7mn ^ = ~^(^ ^1mn~ ^ 

k=0 

k=0 

+ f'( 1 ) _ p(0 ('oil 

fc=0 

i 

_l Vt^^A^ - -G^ 

k=0 

(22) 

where for a generic analytic function f(t) we denned 
as the Z-th coefficient of the Taylor expansion. We now 
show that Eqs. (|2HI22[) constitute a set of recursive re- 
lations to calculate all 7'W and A'^ once all Y and 
A'C ) are known for k < I. Wc first observe that the 
Z-tli derivative of the density matrix p'(t) in t = de- 
pends at most on the (I — 1) derivative of 7' and A' since 
i^/S'CO = [H'(t),p'(t)]. The quantity depends on 
(7', A') implicitly through p'(t) and explicitly through 
the commutator [H'(t), J mn (t)}- Since the Z-th derivative 
of the commutator depends on all (7'W , A'^ ) with k < I 

the quantity F„ is a function of (7'W, A'W) with fc < Z. 
On the contrary, the quantities if', G' depend implicitly 
on (7', A') through p'(t) but they explicitly depend only 
on 7', i.e., there is no explicit dependence on the pair- 
ing potential A'. We therefore conclude that K'"> and 
G'W depend on the 7'W with k < I and on A'W with 
k < I. Finally, from Eq. (flU)) we see that the Z-th deriva- 
tive of the density n' m (i) depends at most on the Z — f 
derivative of 7' and A'. The table below summarizes the 
dependency of the various quantities on the order of the 
derivatives of 7' and A' 





F'V) 


K '{i) 




n '{i) 




k < I 


k < I 


k < I 


k < I 


|A'( fe )} 


k < I 


k < I 


k < I 


k < I 



(23) 



From the above considerations it follows that Eq. (f2"2")) 
with I = can be used to determine A'(°) since the r.h.s. 
depends only on Y' ' = j'(0) and from Eq. (|18| the pref- 
actor [rim — 1] ^ 0. Having A'(°) we can easily calculate 
7'W from Eq. (j2"Tj) with I = since the r.h.s. depends 

only on 7'W and A'<°) and from Eq. (JT7|) ^ 0. 

With 7'^), and A'(°) we can use Eq. (JUJ) with I = 1 
to extract A' (1) , then Eq. ([21]) with Z = 1 to extract 7 /(2) 
and so on and so forth. 



C. Keldysh-Green's function in the Nambu space 

1. Keldysh contour 

Wc now specialize to interacting systems which are ini- 
tially in equilibrium at temperature T = 1//3 and chem- 
ical potential p; such initial configurations are the rele- 
vant ones in quantum transport experiments, see Section 
liTDl ^ From static SCDFTi^ we can choose the initial 
density matrix of the KS system as the thermal density 
matrix of a system described by the equilibrium Hamilto- 
nian ([2]) with H mt — and KS phases 7 and pairing po- 
tentials A, and from the results of the previous section we 
know that such KS system can reproduce the TD bond- 
currents and pairing densities of the interacting system if 
perturbed by TD KS phases 7(f) and pairing potentials 
A(t). Denoting by H s (t) = K(t) + A(t) + At(i) the TD 
Hamiltonian and by p s (t) the TD density matrix of the 
KS system wc then have 

Ps(t) = ^S s (t)e-^»-^St(t) (24) 

where Z = Tr {e - ^(- Hs- ' lJV '} is the partition function 
and S s (t) is the KS evolution operator to be determined 
from i4;S s (t) = H s (t)S s (t) with boundary condition 
S s (0) = 1. The Hamiltonian H s = H s {0) is the equi- 
librium KS Hamiltonian while N is the total number of 
particles operator. It is worth to notice that in general 
[H s , N] 7^ due to the presence of the pairing field. The 
TD expectation value O s (t) of a generic operator 0(t) is 
in the KS system given by^ ~ 36 ' 44 

O s (t) = TY{p s (t)6(t)} = (T K [6{z = t±)}) (25) 

where wc have introduced the short hand notation 



Tr 


m (—if dz H a s (z) 




Tr 







In the above equation 7k is the Keldysh contour— il- 
lustrated in Fig. [1] which is an oriented contour com- 
posed by an upper branch going from to 00, a lower 
branch going from 00 to and a purely imaginary (ther- 
mal) segment going from to —i/3. The operator Tk 
is the contour ordering operator and move operators 
with later contour variable to the left (an extra minus 
sign has to be included for odd permutations of fcrmion 
fields). Finally H^ a (z = t±) — H s (t) where the contour 
points lie on the upper/lower branch at a distance 

t from the origin while for z on the thermal segment 
Hfj,^ s {z = —it) = H s — pN. Thus, the denominator in 
Eq. is simply the partition function Z. In Eq. (|2"5|) 
the variable z on the contour can be taken either on the 
upper it-) or lower (t + ) branch at a distance t from the 
origin. 
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2iT mm t 



z = t. 



0+ 



Z = tj 



-i/3 



FIG. 1: The Keldysh contour 7k described in the main text. 
The contour variable z = t-/t+ denotes a point on the up- 
per/lower branch at a distance t from the origin while z = — it 
denotes a point on the imaginary track at a distance r from 
the origin. In the figure we also illustrate the points 0_ (ear- 
liest point on 7k), 0+ and —i/3 (latest point on 7k). 



2. Keldysh- Narnbu- Green's function 

The KS expectation value O s (t) of an operator 0{t) 
is in general different from the expectation value 0(t) 
produced by the original system. However if 0(t) is the 
KS bond-current operator or the pairing density operator 
the average over the KS system yields exactly the bond- 
current and pairing density of the original system. It 
is therefore convenient to introduce the non-equilibrium 
Nambu-Green's functions (NEGF) from which the ex- 
pectation value of any one-particle operator can be ex- 
tracted. A further reason for us to introduce the NEGF 
is that the equilibrium and time-dependent Bogoliubov- 
deGennes equations can be elegantly derived from them, 
thus illustrating the equivalence between the NEGF and 
the Bogoliubov-deGennes formalisms. The normal and 
anomalous components of the NEGF are defined accord- 
ing trji& 



G 



1 



.„„(z;z') = -(T K {c ma (z)cl a {z')}), 
i 



(27) 



F mn {z;z') = -(T K {c ml (z)c„ t (z')}), (28) 
i 

F m „M = -\{Tk{cI a {z')cI x {z)}), (29) 

where z, z' run on the Keldysh contour 7 K1 34 i 35 i 44 i 47 The 
c operators carry a dependence on the z variable; such de- 
pendence simply specifies their position along the contour 
so to have a well defined action of 7k 4^ The TD bond- 
current and pairing density can be expressed in terms of 
G CT (z; z') and F(z; z') as 

Jmn(fy ^ ^ (7t7X7i,6 ^ ^ ^^ajWmi^ — j ^+) "T" H.C. 



(30) 



3. Equations of motion 



(31) 



The NEGF of the KS system obey the following equa- 
tions of motion 

(4 i " H M (^)| G(z; z') = lS(z - z'), (32) 



G(z;z') 



■X 

dz 



7 l-n fl (z')\=15(z-z'), (33) 



where all underlined quantities are 2x2 matrices in the 
Nambu space with matrix elements l m „ = 1 



and 



<L 



G mn (z; z) 



(*;*') 



.(*';*) 



K„ 



t (z) <5 mn A m (z) 



The matrix elements of H (z) are 

(^±) — T nin G ^ rnn ^ ^ 



A m (t ± ) = A m {t) 
for z = t± on the horizontal branches and 



K 



fi,mn 



( IT) — T mn C ^" 



fJ,S„ 



A m (-ir) = A, 



(34) 



(35) 



(36) 



(37) 



for z = — it on the imaginary track. Since H (— it) is 
independent of r we write H (— ir) = H — /i<r with 
<x m „ = a z l mn and tr z the third Pauli matrix. 

In the next Section we show that the solution of the 
equations of motion is equivalent to first solve the static 
Bogoliubov-deGennes (BdG) equations and then their 
TD version. 



4- Keldysh components and Bogoliubov-deGennes equations 

We introduce the left and right contour evolution ma- 
trices S R ^ L (z) which satisfy 



*^s fl (z) = R,m R (z), 



i^&V) = S L (z')H M (z') 



(38) 
(39) 



with boundary conditions S^ L (0_) = 1. The most gen- 
eral solution of the equations of motion (|32I33[) can then 
be written as 



G(z; z') = S R (z) [9{z\ z')G> + 9(z'; z)G<] S L (z'), 



(40) 
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with G > G < = il and the contour Heavisidc function 
9{z\ z') = 1 if z is later than z' and zero otherwise. Equa- 
tion (l40l) is a solution for all matrices G > 



-il + G\ 

In order to determine G > or G < we use the boundary 
conditions 



G(0_;z') = -G(-ip;z'), 
G>;0_) = -G(z;-ip), 



(41) 
(42) 



which follow directly from the definitions (|27H29[) of the 
NEGF. Using Eq. gD]) one finds G(0_;z') = G^S^-z') 
and G(-iP; z') = S i? (-i/3)G > S i (z') from which we con- 
clude that 



G< = -s^H^cp. 

Similarly, from Eq. (|4"2"j) one finds 

G> = -G<S L (-i/3). 



(43) 



(44) 



Exploiting the fact that H^(— it) = H — [ia_ is con- 
stant along the imaginary track one readily realizes that 
S fl/i H/3) = exp[±/?(H - no)] and hence 



G< = 



l + exp[/3(H -fia)} 



(45) 



From the exact solution (|4T)|) we can extract any observ- 
able quantity at times t > and not only its limiting 
behavior at t — > oo. Below we calculate the different 
components of the NEGF. 

We introduce the eigenstates ^ q , with eigenenergies 
E q . of the matrix H — The vector = 
is a two-dimensional vector in the Nambu space and, by 
definition, satisfies the eigenvalue problem 



E 



T p 1 ^ 



l u q {n) + A m v q (m) = (E q + ^)u q (m), (46) 



E 



Tk m e I7r,m u ? (n) + A* n u 9 (to) 



(£ s - M K(m). (47) 



Due to the presence of the pairing field the components 
u q and v q are coupled and the eigenstates are a 
mixture of one-particle spin- up electron states and spin- 
down hole states. We will refer to the eigenstates <3/ g 
as bogolons. The above equations have the structure 
of the static BdG equations which follow from the BCS 
approximation. 48,49 In our case Eqs. (|46l47p follow from 
SCDFT— and therefore yield the exact equilibrium bond- 
current and pairing density provided that the exact KS 
phases and pairing fields are used. 

Inserting the complete set of eigenstates in Eq. (|4T)]) 
and taking into account Eq. (j45]) we find the following 
expansion for the NEGF 

G(z; *') = < E S R (z)^ q [6(z; z')f>(E q ) 
i 

+9(z';z)f<(E q )]^ q S L (z'), (48) 



where f < {<^>) = 1/[1 + exp(/3w)] is the Fermi function 
and / > (w) = / < (w) — 1. Taking z and z 1 on the real axis 
but on different branches of the Keldysh contour, we can 
extract the lesser and greater component of the NEGF. 
We first notice that for z = t± the contour evolution 
operators reduce to the standard evolution operators, i.e., 
S R (t±) = S(t) and S L (t±) = S}{t) with 



S(0) = 1, 



(49) 



and H(t) = H M±), see Eq. ([35]). Then, in terms of the 
evolved states ^ q (t) = S(t)^ g with components ^ q (t) = 
[u q (t) , v q (t)] we find 



G%i') = G(i T ;4H 



Gf(t;t') 



F^ r (t';t) 



-Gf T (t';t) 



,(50) 



Uq{t) U \{t') U q {t)v\{t') 

v q (t)u\(t>) v q {t)v\{t>) 

where the superscript T in F^ ,T and G^' T denotes the 
transpose of the matrix, see also Eq. ([M]) . The func- 
tions u q (t) and v q (t) can be determined by solving a cou- 
pled system of first-order differential equations. From 
Eq. (gl]) it follows that 

i^-u q (m, t) = T mn e tlmni - t] u q (n, t) + A m (t)v q (m, t), 



(51) 



i—v q (m,t) = - ^T„ m e l7 " m(t) -u g (n,t) + A* m (t)u q (m,t), 

n 

(52) 

which have the structure of the TD BdG cquationsi 26 i 50 
As in the static case, however, the solution of Eqs. ([SH - 
IS!?]) yields the exact densities and not their BCS approx- 
imation. 

We notice that for the KS system to reproduce the 
time-independent densities of an interacting system in 
equilibrium it must be 



A m (t)=e- 2 ^A r , 



(53) 



for which one finds the solutions u q (t) = e ^i+^Uq 
and v q {t) = e~ l ^ Eq ~ IJ ^ t Vq. The above time-dependence of 
the pairing field is the same as in the BCS approximation. 

Using Eq. ([50]) the retarded (R) and advanced (A) 
NEGF are 

G R/A {t;t') = ±6(±tTt')[G > (t;t')-G < (t;t')] 

= T i8(±tTt')S{t)S}(t'), (54) 

with components 



G R m / n A (t;t') = 



G 
F 



R/A (t- 1') 

R/A 



(f;t) 



(*;f) 



xpA/R 
-i- nm 

A/R (t'-t) 



G 



(55) 



8 



It follows that G^(t; t') can also be written as 

G${t;t') = G R (i;0)G^(0;0)G A (0;i')- (56) 

D. Application to quantum transport 

We here apply the above formalism to systems de- 
scribed by a = l,...,Af bulk superconducting leads in 
contact with a central region C which can be, e.g., a 
quantum dot, a molecule or a nanostructure. Assuming 
no direct coupling between the leads the Hamiltonian 
H „ is written in terms of its projections on different sub- 
spaces as 

AT 

—fi,aa 



H, 



H 



p,,CG 



AT 



H,, C J, (57) 



a=l 



where H aa describes the a-th lead, H^ cc the nanos- 



tructure C and H aC 



the coupling between 
lead a and C. We assume region C to be a constric- 
tion so small that the bulk equilibrium of the leads is 
not altered by the coupling to C . Furthermore we con- 
sider time-dependent perturbations which correspond to 
the switching on of a longitudinal electric field in lead 
a. The time to screen the external electric field in the 
leads is in the plasmon time-scale region. If we are in- 
terested in external fields which vary on a much longer 
time-scale it is reasonable to expect that the leads remain 
in local equilibrium. Therefore the coarse-grained time 
evolution of the system can be described by the following 
±) = H(i) 



TD Hamiltonian (t 



H aa (<) = exp (-ifxta z ) H aQ (0) exp (i/j,ta z ) . 



H aC (t) = exp [i / diU a (t)a z H aC (0) 



(58) 



(59) 



Hoa(*) = HcJi)] 1 - (60) 

We do not specify the time dependence of H cc (i) since 
it can be any, see below. The TD field U a (t) is the sum of 
the external and Hartree field and is homogeneous, i.e., 
it does not carry any dependence on the internal struc- 
ture of the leads, in accordance with the above discus- 
sion. It has been shown that for macroscopic leads the 
assumption of homogeneity is verified with rather high 
accuracy^ 

As for the case of normal leads the equations of motion 
for the Keldysh-Green's function can be solved by an em- 
bedding procedure. We define the uncontacted Green's 
function g which obeys the equations of motion (|32I33[) 
with H„ aC = H„ Ca = and the same boundary con- 
ditions as G. Then, the equation of motion for G cc 
projected onto regions CC takes the form 

r.t 

dz ' 



— Xcc ~ H M cc (z) >G cc (z; z') = l cc 6(z 



vj dzUz;z)G cc (z; z'),(61) 



where the embedding self-energy is expressed in terms of 
g as 

N 

S(z;z') = £S Q (z;z') 

a=l 
Af 

= Y,^,c a ^)S a Jz;z')U^ aC (z'). (62) 



The above equation of motion is defined on the Kcldysh 
contour of Fig. [T] Converting Eq. (|6I|) in equations for 
real times results in a set of coupled equations known as 
Kadanoff-Baym equation a 34 ' 52 — recently implemented 
to study transient responses of interacting electrons in 
model molecular junctionsi 51 i 57 The use of the Kadanoff- 
Baym equations to address transient and relaxation ef- 
fects in other contexts has been pioneered by Schafer, 58 
Bonitz et al.^ and Binder et al.i^S 

The importance of using an uncontacted Green's func- 
tion g with boundary conditions (|41I42[) for a proper de- 
scription of G > (t\ t') at finite times has been discussed 
elsewhere in the context of transient regime a 36 ! 51 and it 
has been shown that it leads to coupled equations be- 
tween the Keldysh-Green's function with two real times 
and those with one real and one imaginary time. 

In the next Section we propose a wave-function based 
propagation scheme to solve Eq. (I5T1) for TD Hamiltoni- 
ans of the form (|58HG0D . 



III. NUMERICAL ALGORITHM 

We consider semi-infinite periodic leads with a super- 
ccll of dimension N^ ell for lead a. The projected Hamil- 



tonian H 



0,aa 



H QQ (0) can then be organized as follows 



H 



■O.aa. 



ha 


£o 


% 




h a 


to 


o„ 


tl 


h a 



(63) 



where h a is the 2N" ell x 2N" ell Nambu Hamiltonian of 
the supcrcell with matrix structure 



h, 



A* 



(64) 



while t a describes the contact between two nearest neigh- 
bor supercells. Since the pairing field is local the off- 
diagonal terms of t a are zero and therefore the general 
structure of the hopping matrix is 



ta Q 

o„ -tl 



(65) 



The matrices e Q , A Q and t a in and t a have the di- 
mension of the unit cell, i.e., N^ oll x N" ea . In particular 
A Q is a diagonal matrix. 
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A. Calculation of initial states 

Given the above structure of the leads Hamiltonian 
the eigenstates of H — fig_ can be grouped in scattering 
states with incoming bogolons from lead a = 1 , . . . , Af 
and Andreev bound states (ABS). 



1. Scattering states 

The lead a is characterized by energy bands E™(p) 
with v = 1 7 . . . ,2N" ell and p <E (0, n). For a given p 

I 



the energies E"(p) are the solutions of the eigenvalue 
problem 



lp 



■tie-* 



fJ-°Lc 



jja 



K(P)U" (66) 



with U" p the Nambu-Bloch eigenvectors. We write the 
index of the localized orbital ip m as m = s,j, a; here s 
labels the orbital within the supercell, j the supercell and 
a the lead. The index s runs between 1 and N" e ^ while 
the supercell index j = 0, . . . , oo. The scattering state for 
an incoming bogolon from lead a has the general form 



( U? p (s)e-™ + £ R% tp W™ p (s) e 



< P (™) = { 



m = s,],a 
meC 

m = s,j,(3 ^ a 



(67) 



with reflection coefficients R and transmission coefficients 
T. The momenta q"fi p (for all leads /3 including f3 = a) 
are associated to states with energy E = E"(p) and can 
therefore be obtained from the roots of 



Det[/ig + tpe 



tie 



El, 



0. 



(68) 



The above equation admits, in general, complex solu- 
tions for q. In Eq. (|67[) the sums over p run over 
real solutions q for which the sign of the Fermi velocity 
v p (l) = dEjj (?) I i s opposite to the sign of the Fermi 
velocity {p) of the incoming bogolon and over all com- 
plex solutions q for which Im[g] > (evanescent states). 
Once the q%ff p arc known the Bloch state is sim- 

ply the eigenvector with zero eigenvalue of the matrix 

3 



h 



toe 



"I,. 



tje-^.p 



Elg. For the cal- 



culation of the reflection and transmission coefficients as 
well as of the amplitude c (to) in the central region 
we extended a recently proposed wave-guide approach^ 
The method is based on projecting the Schrodinger equa- 
tion (H — /i<x)<I' = E^> onto the central region and onto 
all the supercells in contact with the central region, i.e., 
with j = 0. The projection onto a j = supercell leads 
to an equation which couples the amplitude of V E' in j = 
with that in j = 1. Exploiting the analytic form of the 
eigenstate in Eq. (|57|) the amplitude in the leads can 
entirely be expressed in terms of the unknown i?'s and 
T's for all j. In this way the equations can be closed 
and the problem is mapped into a simple linear system 
of equations for the unknown R"P p , ' and c( m )- 



2. Andreev bound states 

The presence of a gap in the spectrum of the supercon- 
ducting leads may lead to the formation of localized ABS 
within the gap. The procedure to calculate the ABS is 
slightly different from the one previously presented since 
the ABS energy is not an input parameter and the ABS 
state is normalized to 1 over the whole system. The en- 
ergy Eb of an ABS VE^ is outside the lead continua. Pro- 
jecting the Schrodinger equation (H — /-ig^^b = E^b 
onto different regions and solving for the projection ^b,c 

in region C one finds (H.q cc (Eb) — fia cc )^b,c = Eb^b,c 
where 



H 



o,cc 



(E)=R 



o,cc 



1 



0,Ca 



^-(H _ QQ -/iCT QQ )" 



-H 



■0,qC- 



(69) 



The ABS energies Eb can then be extracted from the 



— El cc ] = and the 



roots of Bct[H^ cc (E) 

eigenvector with zero eigenvalue of II.Q S cc (Eb) — p(L cc 
EbXcc i s proportional to the projection *f?b,c of the ABS 
in region C . We call C& the unknown constant of propor- 
tionality. As for the scattering states we can construct 
the ABS everywhere in the system according to 



*b(m) 



p 



(70) 



m e C 



The momenta g" p and Bloch states W b a are calculated 
in the same way as for the scattering states. By definition 
all momenta have a finite imaginary part and the sum in 
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Eq. ()70|) runs over those with a positive imaginary part. 
The constants B^ p can be simply obtained by project- 
ing the Schrodinger equation (H — /la^b = onto 
the supercells in contact with region C , i.e., with j = 0. 
The resulting equation couples the amplitude of ^b in 
j = with that in j = 1 and with the known amplitude 
Cb^b.c( m )- Exploiting the analytic form of \I/f, in the 
leads the amplitude in j = 1 can entirely be expressed 
in terms of the constants CbB" p thus yielding a linear 
system of equations for each lead. Once the CbB£ p are 
known the constant of proportionality Cb is fixed by im- 
posing that the ABS is normalized to 1. This can be 
easily done since the sums over j are geometrical series. 



perfect transparent boundary conditions at the interfaces 
between region C and leads a. Projecting Eq. (|75p onto 
lead a and iterating one finds 



$(m+i) = „m+U(o)_ 



iSH, 



~ T (m-j) 
C 



3=0 



X $ 



(m+l-j) 
C 



where we have defined the propagator 



C 



,(76) 



(77) 



Embedded Crank-Nicholson propagation 
scheme 



To propagate the generic eigenstate of H — /i<x we 
extend the embedded Crank-Nicholso n 37,38 scheme to su- 
perconducting leads. The equations of motion (|51I52[) 
can be written in a compact form as 



i— *(f) =H(f)#(t), *(0) = * 



(71) 



where the components of the TD Hamiltonian are given 
in Eqs. (|58H60[) . We first perform the gauge transfor- 
mation ty a (t) = cxp[— i^CT QQ t]$ Q (t) for the projection of 
the state onto lead a and ^c(t) = &c{t) for region C. 
The state $(£) obeys the equation 



ij t Ht) = H(t)*(f), *(0)=* (72) 



with 



Haa W - Haa(°) _ M^aa 



(73) 



and made use of the fact that H QQ (i) = H QQ is time- 
independent. The time-dependence of the contacting 
Hamiltonian can be easily extracted from Eq. (|74[) and 
reads 



cxp 



(■ (m+l) 



where we have defined 



(m) 

CXp l/i a 'CT^ 



dtU a (t). 



H QC (0), 
(78) 

(79) 



At this point comes a crucial observation which allows 
for extending the propagation scheme of Rcfs. I37H381 to 
the superconducting case. Since the pairing field is local 
in the chosen basis the off-diagonal part of the contacting 
Hamiltonian is zero and hence H Ca cr aa = <r cc H Ca . It 
follows that Eq. ((78)) can also be rewritten as 



fi C*0 



. aC ~ H aC (0)- 



(m+l) 

exp I g^cc 



(m) 

exp l 'g^jc 



Hac(*) = ex P 



lit 



dtU a (t) 



H ttC (0) 



and H cc (t) 



(74) 

The advantage of the gauge 
transformed equations is that the lead Hamiltonian is 
now independent of time. We discretize the time as 

t rn = 2mS and define $ (m) 
l 

2 



H(i 



m+l, 



$(t m ) and H (m) = 
The differential operator in Eq. 



(I72[) is then approximated by the Cayley propagator 



H aC (o)iL m ), 



(80) 



which implicitly define the matrices z_a = (z_a )*• 
Next we project Eq. (|75p onto region C and use Eq. 
(|76l) to express the $ Q at a given time step in terms of 
the $c a t a ll previous time steps. The resulting equation 
is 



Ice 



Icc - tffiS^ 



( s * n) + Mi a m) ) (81 



1 + z5H (m) ) = ( 1 - zffl 1 ^ )$(»»). (75) 



r x ( m ) 



The above propagation scheme is known as Crank- 
Nicholson algorithm and it is norm-conserving and ac- 
curate up to second order in 5. As the matrix H is in- 
finite dimensional the direct implementation of Eq. (|75[) 
is not possible. A significant progress can be done using 
an embedding procedure which, as we shall see, entails 



and contains only quantities with the dimension of region 
C. We emphasize that Eq. (|81l) is an exact reformulation 
of the original Eq. (|75|) but it has the advantage of being 
implcmentable. Indeed, exploiting the result in Eq. 



the boundary term and memory term M a read 



Sz^U Ca (0) g™ (l aa + k) $i 0) , (82) 
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x ( $p J ' + <i>p J ' 



while the effective Hamiltonian is given by 



T ~ T (ni) T ~ T ("T.) 

H cff = H cc - 



(83) 



(84) 



where the embedding matrices Q*-™-' have twice the di- 

— ^ a 

mension of region C and are defined according to 

Qi m) =H Ca (0)^ _Z_H ac(0 ). (85) 

In Appendix 1X1 we describe a recursive scheme to calcu- 
late the embedding matrices. In Appendix IB1 wc further 
show that the boundary term Sa can be expressed in 
terms of the Q 's thus rendering Eq. (|8ip a well defined 
equation for time propagations. 

In the next Section we apply the numerical scheme to 
UF-JNJ model systems and obtain results for the TD 
densities and currents. 



describes the lead a = L,R, 

H LC (t)=t LC e^o^u L (t')^2cl aL c 0r7 + ll.c. (87) 



H RC (t) = tnce'fo rft 'tM*') £ cl R c Na + H.c. 



accounts for the coupling between region C and the leads, 
and 



JV-l 



N 



H C c =tcYl E( £ 

+H.c.)+ ec E E 



m— a 



m— cr 



(89) 

is the Hamiltonian of the chain with N + 1 atomic sites. 
The currents J L (t) = J l,o(0 and J R (t) = Jjv,ofl(0 
through the bonds connecting the chain to the left and 
right leads arc obtained from Eq. (|3U1) and Eq. ([5U)) and 
read 



-itLce 



ilhc(t) 



Y,f>(E g )v q (Q,t)v* q (OL,t) 



H.c. 



(90) 



IV. REAL-TIME SIMULATIONS OF S-D-S 
JUNCTIONS 

Due to the vast phenomenology of S-D-S junctions it 
is not possible to address these systems in a single work. 
Furthermore the analysis of the time-dependent regime 
is generally more complex than that in the Josephson 
regime and it is therefore advisable to first gain some 
insight by investigating simple cases. Our intention in 
this Section is to demonstrate the feasibility of the prop- 
agation scheme and to present genuine TD properties of 
simple model systems. 

We consider a tight-binding chain (region C) with 
nearest neighbor hopping tc and on-site energy ec con- 
nected to a left (L) and right (R) wide-band leads. The 
a = L, R lead is described by a semi-infinite tight-binding 
chain with nearest neighbor hopping t a and a constant 
pairing field A Q , and is coupled to the a end-point of 
the central chain through its surface site with a hop- 
ping tea = t a c- The system is initially in equilibrium 
at temperature T = and chemical potential fj, = 
and driven out of equilibrium by a TD bias voltage 
U a (t) applied to lead a at positive times. From Sec- 
tion III D[ the Hamiltonian for this kind of systems read 
H(t) = E a (H aa (t) + H aC (t) + H Ca {t)) + H CC where 

OO 

H aa (t) = t a J2J2(c] +1<ja c ]aa +ll.c.) 

j=0 a 

+ (e-^'A^t 4+H.c.) (86) 



Mt) = -HI 



RC 1 



-ilRc{t) 



J2.f < (E q )u q (0,t)u* q (0R,t) 



L q 



q 



f>(E q )v q (0R,t)v*J0,t) 



H.c. 



(91) 



where J a c(t) — i f£ dt'U a (t') and the sum over q runs 
over all ABS and scattering states. Similarly, the pairing 
density P m (t) on an arbitrary site of the chain is obtained 
from Eq. (J3T]) and Eq. (J5DJ and reads 



Pratt) = J^/^KKtKKt) 

q 



,2ie c t 



(92) 



We will write the pairing field as A Q = £ Q e lXa A and 
measure energies in units of A, times in units of h/A 
and currents in units of \e\A/h, with |e| the absolute 
charge of the carriers. Since we consider wide-band leads 
with t a 3> t a c,tc and the chemical potential is set to 
zero the results depend only on the ratio T a = 2t^ c /t a 
(tunneling rate) and not on t a c and t a separately. In 
the following we therefore specify the value of T a only. 
In practical calculations the longitudinal vector p G (0, tt) 
of the scattering states, see Eq. (joTf . is discrctized with 
N p mesh points and only states with energy within the 
range (fj, — A, [i + A) are propagated in time. We will call 
Np a the number of scattering states from lead a that 
are propagated. The cutoff A is chosen about an order 
of magnitude larger than the typical energy scales of the 
problem, i.e., U a , T a , A Q , t c , ec- 
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A. The single-level quantum dot model 

The single-level quantum dot (QD) model corresponds 
to a central chain with only one atomic site (N = 0). For 
A L = A R = (N-QD-N) the TD response of this system 
has been investigated by several authors and an analytic 
formula for the TD current is also availablei 36 i 62 ' 63 Scarce 
attention, however, has been devoted to the system with 
one superconducting lead2S (N-QD-S) and to the best 
of our knowledge the only available results when both 
leads are superconducting (S-QD-S) have been published 
in Ref. HI 



N-QD-S model under DC 



We first consider the N-QD-S case schematically illus- 
trated in Fig. [21a). To highlight the different scattering 
mechanisms we shift the central level by ec = 0.5, choose 
weak couplings to the leads Fl = Fr = 0.2, and drive 
the system out of equilibrium by applying four different 
biases Ul = 0.3, 0.6, 0.9, 1.2 to the left normal lead. For 
biases in the subgap region, i.e., Ul < A# = 1, transport 
is dominated by Andreev reflections (AR). In Fig. H£b) 
we show the currents J Lit) and Jr{£) of Eqs. (|9QI9ip . 
For Ul =0.3 < ec the AR are strongly suppressed since 
electrons at the left electrochemical potential \Il = Ul 
have just enough energy to enter the resonant window 
(e c - 2r, e c + 2r), where 2r = T L + T R . Resonant AR 
can occur for Ul > ec and constitute the dominant mech- 
anism for electron tunneling. This is clearly visible in the 
second panel of Fig. HJb) where the steady-state values 
of Jr for Ul = 0.6 and Ul = 0.9 are approximatively 
the same. At larger biases Ul = 1.2 > Ar electrons can 
also tunnel via standard quasi-particle scattering and the 
steady-state current increases. This interpretation is con- 
firmed by the behavior of the pairing density Poit) on the 
QD, third panel of Fig. HJb). For times up to ~ 5 the 
pairing density decreases since pre-existent Cooper pairs 
in lead R move away from the QD. However, while |Po(t)| 
remains below its equilibrium value at Ul = 0.3, for all 
other biases, Ul > ec, |-Po(*)| increases after t ~ 5, mean- 
ing that a Cooper pair is forming at the interface. We 
also notice that the values of \Po(t —> oo)| for Ul = 0.9 
and Ul = 1.2 are very close while the corresponding cur- 
rents Jr differ appreciably. This is again in agreement 
with the fact that electrons with energy larger than Ar 
do not undergo AR and thus no extra Cooper pairs are 
formed. Finally we observe that the transient regime is 
longer in the N-QD-S case than in the N-QD-N case, see 
inset in panel 2 and 3 of Fig. HJb), as also pointed out in 
Ref. HI 



2. S-QD-S model under DC bias 

We now turn to the more interesting case in which 
the QD is connected to a left and right superconduct- 
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FIG. 2: (a) Schematic of the transport set up. A single 
level QD with on-site energy ec = 0.5 is weakly connected 
(r_L = r_R = 0.2) to a left normal lead and a right supercon- 
ducting lead. In equilibrium both temperature T and chemi- 
cal potential /i are zero. The system is driven out of equilib- 
rium by a step-like voltage bias Ul — 0.3, 0.6, 0.9, 1.2 in the 
normal lead. For Ul < A_r the dominant scattering mecha- 
nism is the AR in which an electron is reflected as a hole and 
a Cooper pair is formed in lead 7?. (b) Time-dependent cur- 
rent at the left interface (first panel), right interface (second 
panel) and absolute value of the pairing density on the QD 
(third panel). The insets show the TD current for the same 
parameters but A_r = 0, i.e., for a normal R lead. The results 
are obtained with a time-step 8 — 0.05, cutoff A = 6 and a 
number of scattering states N Pi l — 1070, N p< r = 1056. 



ing lead (S-QD-S), sec Fig. 
metric couplings Tl = Fr 
fields Al = A R e ix = ^ 
but different phase. 



G3a). We focus on sym- 
= r = 1 and on pairing 
e" A with the same magnitude 
This system always support two 



Andreev bound states (ABS) in the gap. Their en- 
ergy can be obtained analytically from the solution of 
Det[Ho ff cc (£') - /xctcc ~ E hjc\ = ( see Section liTTX2|) 
which, in terms of the dimensionless variables x = E/ A, 
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FIG. 4: (a) Discrete Fourier transform of Jl(£) in arbitrary 
units [the curves corresponding to bias Ul — O.n are shifted 
upward by 0.7(n — 3) while that corresponding to bias Ul = 
1.0 is shifted upward by 2.8. (b) Values of the average current 
for biases in the subgap region, (c) ABS contribution to the 
current Jzjf) for biases Ul = 0.2, 0.3, 0.4, 0.5, 0.6 [the curves 
corresponding to bias Ul = O.n are shifted upward by 0.8(71 — 
2)]. The numerical parameters are the same as in Fig. 



FIG. 3: (a) Schematic of the S-QD-S model with F L =r R = 
1.0, Al = Ar = 1, and ec — 0. This system admits two ABS 
in the gap. The ABS energy depends on the superconducting 
phase difference \ as illustrated in the inset, (b-c) Time- 
dependent current JL{t) at the left interface as a function of 
time for (b) Ul = 3.0, 2.0, 1.0 [the curves corresponding to 
bias Ul — n.O are shifted upward by 0.3(n— 1)] and (c) Ul = 
0.5, 0.4, 0.3, 0.2 [the curves corresponding to bias Ul = O.n 
are shifted upward by 0.6(n — 2)]. The results are obtained 
with a time-step S = 0.05, cutoff A = 12.1, and a number 
of scattering states N Pi l = N p> r = 768 for panel (b) and 
6 = 0.05, A = 4, N P , L = N p ,r = 788 for panel (c). 



7 = I 1 / A and e = (ec — m)/A, reads 

„2n , 7 



r 



2 «V 



1-x 2 



(93) 



where a = y 1+c 2 osx and varies in the range (0, 1). In Fig. 
[3]Ja) we plot the solutions of Eq. (|9"3"| as a function of x 
for ec = H = 0. In equilibrium and at zero temperature 
one ABS is fully occupied and the other is empty. At time 
t = a constant bias Ul is applied to the left lead. In 
Fig. [2Kb) we display the TD current at the left interface 
Ji(t) for x = an d Ul = 3, 2, 1. After a transient the 
current oscillates in time with period Tj = 2tt/(2Ul), 
as expected. For Ul > 2 the S-QD-S system behaves 
similarly to a macroscopic Josephson junction with an 
almost pure monochromatic response, albeit the average 
value Jdc of the current over a period is different from 
zero. For Ul = 1 < 2A, i.e., in the subgap region, the 
transient regime becomes much longer and Jl (t) deviates 



from a perfect monochromatic function. At Ul = 1 the 
dominant scattering mechanism is the single AR. 

As discussed in Ref. [TH the presence of the resonant 
level modifies substantially the Jdc — V (V — Ul — Ur) 
characteristic and for T = 1 the subharmonic gap struc- 
ture is almost entirely washed out. However, a very rich 
structure is observed in the TD current. In Fig. (3Jc) 
we display J L (t) for biases U L = 0.5, 0.4, 0.3, 0.2. The 
charge carriers undergo multiple AR (MAR) before ac- 
quiring enough energy and escaping from the QD. The 
dwelling time increases with decreasing bias and the tran- 
sient current has a highly non-trivial behavior before the 
Josephson regime sets in. From the simulations in Fig. 
[3jc) at bias Ul = 0.2 the propagation time t = 250 is 
not sufficient for the development of the Josephson oscil- 
lations. We also observe that the smaller is the bias the 
larger is the contribution of high-order harmonics, which 
is in contrast with one would naively expect from linear 
response theory. 

In Fig. HJa) we display the Fourier transform of 
Ji(t) — Jdc in the Josephson regime. Replica of the main 
Josephson frequency wj = 2LJ are clearly visible for 
Ul < A. The values of Jd c as obtained from time propa- 
gation are reported in Fig. Ub) and are consistent with 
a smeared sub-harmonic gap structure. 

From the curves Jl (t) it is not evident how to estimate 
the duration of the transient time. We found useful to 
look at the contribution of the ABS, Jl,abs> to the total 
current Jl, since Jl,abs(£ — > °°) = 0. This quantity is 
evaluated from Eq. (j90|) by restricting the sum over q to 
the ABS and is shown in Fig. H£c). ABS play a crucial 
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FIG. 5: Time-dependent current at the right interface Jr 
(first panel) as well as the density no (second panel) and pair- 
ing density |Po| (third panel) on the QD. The curves from bot- 
tom to top corresponds to a switch-off time = 5tt + nir/8, 
with n = 0, 1, 2, 3, 4. Since the bias is JJh = 1 the ac- 
cumulated phase difference \ a t the end of the pulse is 
yXn) _ 2v£g — mr/4. For the switch-off time t^g the curves 
of Jr are shifted upward by 0.3n, those of no by 0.5n and 
those of | Po| by 0.2n. The results are obtained with a time- 
step 5 = 0.05, cutoff A = 12.1, and a number of scattering 
states N P , L = N P>R = 768. 



role in the relaxation mechanism as we shall see in the 
next Section. 



3. S-QD-S model under DC pulses 

As mentioned in the introduction the possibility of em- 
ploying UF-JNJ in future electronics rely on our under- 
standing of their TD properties. In the previous Section 
we studied the transient behavior of a S-QD-S system 
under the sudden switch-on of an applied bias. Equally 
important is to study how the system responds when the 
bias is switched off. We therefore consider the same S- 
QD-S model as before with Yl, = Tr = 1, ec = 0, 
Al = Ar = 1 initially in equilibrium at zero temper- 
ature and chemical potential. At time t = a constant 
bias Ul = 1 is applied to lead L until the time t a g at 
which the bias is switched off. How does the system re- 
lax? In Fig. [5] we show the current Jr at the right 
interface as well as the density uq and pairing density 
| Pol on the QD for switch-off times t^g = 5n + rnr/8 
with n = 0, 1, 2, 3, 4. Despite the fact that the switch- 



off times are all very close [t^ ~ 15.71 and ~ 17.28] 
the system reacts in different ways and actually relaxes 
only in one case. The strong dependence on t a s is due 
to the two ABS in the gap. Similarly to what happens 
in normal systems^ the asymptotic (t — > oo) form of the 
density on the QD is 

n (t) - no.cont ~ Yl & C0S (( £ ABS - e ABs)*)' ( 94 ) 

ij 

where e^ BS , i = 1,2, are the ABS eigenenergies of the 
Hamiltonian after the bias has been switched off and 
"o,cont is the contribution of the continuum states to 
the density. The coefficients — fji are matrix ele- 
ments of the Fermi function f(H(0)) calculated at the 
equilibrium Hamiltonian and depend on the history of 
the applied biasi 65 i 66 Contrary to the normal case, how- 
ever, the energy of the ABS depends on when the bias 
is switched off since after a time i ff the phase difference 
X changes from zero to 2[/z,£ ff- This fact together with 
Eq. (|M|) explains the persistent oscillations at different 

frequencies. Indeed \^ = 2[/£,i^ = n7r/4 and from 

Fig. Ha) we see that [e^ s ( X (n) ) - 4^ S (x (n) )] varies 
from ~ 1.08 to zero when n varies from zero to 4. The 
amplitude of the oscillations as well as the average value 
of the density no, however, do not depend only on \ but 
also on the history of the applied bias. Two different 
biases Ul{P) and U' L (t) yielding the same phase differ- 
ence X = 2 J to!! dTU L (r) = 2 J * otf <ItU' l {t) give rise to 
different persistent oscillations, albeit with the same fre- 
quency. 

From the results of this Section we conclude that for 
devices coupled to superconducting leads a small differ- 
ence in the switch-off time of the bias can cause a large 
difference in the relaxation time of the device. This prop- 
erty may be exploited to generate zero bias ac currents 
of tunable frequency. 



4- S-QD-S model under AC bias 

The time-propagation approach has the merit of not 
being limited to step-like biases as it can deal with 
any TD bias at the same computational cost. Of spe- 
cial importance is the case of ac biases where a mi- 
crowave radiation U r sin(w r i) is superimposed to a dc sig- 
nal V = Ul — {Jr. The study of UF-JNJ in the presence 
of microwave radiation started with the work of Cuevas 
et ali^i who predicted the occurrence of subharmonic 
Shapiro spikes in the Jd c — V characteristic of supercon- 
ducting point contacts. Later on Zhu ct ali£& extended 
the analysis to the S-QD-S model and discuss how the 
ABS modify the Jd c — V characteristic. The replicas of 
the Shapiro spikes have been experimentally observed^ 
and can be explained in terms of photon-assisted mul- 
tiple Andrcev reflections. Using a generalized Floquct 
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FIG. 6: (a) TD current at the left interface for Ul = 0, 
I7 r = 0.05ct; r with oj r = 0.5, 1.08 [the curve is shifted upward 
by 0.4], and 1.5 [the curve is shifted upward by 0.8]. (b) 
ABS and continuum contribution to the total current in the 
resonant case w r = 1.08, U r = 0.05oj r and Ul = 0. (c) Pairing 
potential on the QD for the same parameters as in panel (b). 
The results are obtained with a time-step S — 0.05, cut-off 
A = 4, and a number of scattering states N Pi l = N p ,r = 788. 



formalism one can show that in the long-time limit^ 7 - 

J L (t) = £ JZ(V, 7 , u T y^ +"^ 4 (95) 

mn 

where 7 = U T /u> T and Uj = TV is the Josephson fre- 
quency. The calculation of is, in general, rather com- 
plicated and to the best of our knowledge the full TD 
profile of J Lit) as well as the duration of the transient 
time before the photon-assisted Josephson regime sets in 
have not been addressed before. 

We here consider the S-QD-S model with Tl = = 1, 
£c = 0, Al = A 7? ; = 1 under a dc bias and in the presence 
of a superimposed microwave radiation f/j,(t) = Ul + 
U T sin(w r i) and Ur = 0. In Fig. EJa) we display the TD 
current at the left interface for fixed 7 = U r /oj r = 0.05 
and different values of the frequency u; r = 0.5, 1.08, 1.5. 
The first striking feature is the occurrence of a transient 

resonant effect at uj r = 1.08 ~ uabs = £ abs — c abs- At 
the resonant frequency the amplitude of the oscillations 
increases linearly in time till a maximum value ~ 0.3. 
The Fourier decomposition (not shown) reveals that the 
peak at to = 1.08 splits into two peaks, one above and one 
below 1.08, which is consistent with the observed beat- 
ing. The effect is absent at larger (w r = 1.5) and smaller 



(w r = 0.5) frequencies for which the amplitude of the 
oscillations remains below 0.05 and two main harmonics, 
one at w r and the other at wabs > are visible in the Fourier 
decomposition (not shown) . The peak at to = wabs is due 
to a transient excitation with a long life-time and cannot 
be described using Floquet based approaches. 

The ABS play a crucial role in determining the TD 
profile of Jl at the resonant frequency. The total current 
J Lit) = JL,cont(t) + Jl, abs it) is the sum of the current 
^L.cont coming from the evolution of the continuum states 
and the ABS current Jl. ABsit)- These two currents are 
shown in Fig. [f^b) from which it is evident that ABS 
carry an important amount of current not only in the dc 
Josephson effec t 30 ! 70 but also in the transient regime. In 
Fig. [6^c) we show the pairing density on the QD for the 
resonant frequency uj r = 1.08. 

In the presence of an external bias the ABS contribute 
to the current only in the transient regime. The duration 
of the transient is investigated in Fig. [7] where we show 
Jr,abs for dc biases with a superimposed microwave 
radiation described by J7z,(t) = Ul + U T sin(o; r t) , with 
U t = 0.05w r , u T = 1.08, and U L = 0.0, 0.03, 0.1, 0.3. 
The interplay between the ac Josephson effect and the 
resonant microwave driving leads to complicated TD pat- 
terns for small Ul- Increasing Ul the life-time of the 
quasi ABS decreases resulting in a fast damping of the 
oscillations, see Fig. [7] with Ul = 0.3. 



B. Long atomic chains 

We consider a chain of N + 1 = 21 atomic sites with 
onsite energy ec = and nearest neighbor hopping tc = 
1, see Eq. (|89]) . symmetrically coupled, Tl = Fr = 
r, to superconducting electrodes with |A^| = |Aij| = 
A. In the limit of long chains one can prove that the 
current phase relation (at zero bias) is linear if tc = 
r/2i 30 ' 70 This is the so called Ishii's sawtooth behavior— 
and is due to perfect AR. To better visualize the MAR 
in the transient regime we therefore choose tc = T/2. In 
equilibrium there are 16 ABS in the gap. At time t = 
the system is driven out of equilibrium by a dc bias Ul 
applied to lead L. 

In Fig. [S] we display the contour plot of the cur- 
rents J n ,n+i(t) along the bond (n,n + 1) of region C 
as a function of time for different values of Ul = 
2A/4, 2A/3, 2A/2. The MAR pattern is illustrated 
with black arrows. There is a clear-cut transient sce- 
nario during which electrons undergo n AR before the 
ac Josephson regime sets in, with n = Ul/2A. At every 
AR the current increases since the electrons are mainly 
reflected as holes and holes as electrons. The same nu- 
merical simulation in a normal system would have given a 
current in region 1AR smaller than the current in region 
OAR. 

For the same system parameters we also considered a 
dc bias Ul = 0.8 for which the dominant scattering mech- 
anism is the 3-rd order AR. The contour plot of the bond 
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FIG. 7: ABS contribution to the current at the right interface for dc biases with a superimposed microwave radiation described 
by U L (t) = U L + U, sm(w r t), with U r = 0.05^ r , uj r = 1.08 and U L = 0.0, 0.03, 0.1, 0.3. The system is the same as in Fig. |g] 
with At = Ar = 1, r L = r R = 1 and e c = 0. The time-step is S = 0.05. 
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FIG. 8: TD picture of MAR. A chain of 21 atomic sites is symmetrically connected with Tl — Fr = 2tc = 2 to two identical 
superconducting leads with Al = Ar = 1. A dc bias Ul = 2A/n, n = 4, 3, 2, is applied to lead L at time t = 0. The panels 
show the contour plots of the bond-current J n ,n+i(t) across the atomic bonds of region C. The results are obtained with a 
time-step S = 0.05, cut-off A = 4 and a number of scattering states N Pi l = N Pt R = 1232. 



current is displayed in the top- left panel of Fig. [S] and is 
similar to the case Ul = 2A/3 of Fig. [5] A new scat- 
tering channel docs, however, open if a microwave radia- 
tion of appropriate frequency is superimposed to Ul ■ We 
therefore applied an ac bias Un(t) = U T sin(o; r t) to lead 
R and choose w r to fulfill 2Ul + oj r = 2A, i.e., uj r = 0.4. 
In Fig. [9] we report the contour plot of the bond-current 
for different values of U x = 0.0, 0.1, 0.3, 0.5. At U T ^ 
the right-going wave-front reduces its intensity just af- 
ter crossing the bond 10 due to scattering against the 
left-going wave-front from lead R, see the characteris- 
tic A-shape in the bottom-right panel. When the right- 
going wave-front hits the right interface the bond cur- 
rent sharply increases. Furthermore, the larger is U r the 
shorter is the transient regime. This can be explained as 
follows. At large U r the dominant scattering mechanism 
is the one in which an electron from lead L and energy 
Ul is reflected as a hole and at the same time absorbs 
a photon of energy u> T . The energy of the reflected hole 
is 2Ul + u>i = 2 A, no extra AR are needed for charge 
transfer and the photon-assisted Josephson regime sets 
in. 



V. CONCLUSIONS AND OUTLOOKS 



In this paper we proposed a one-particle framework 
and a propagation scheme to study the TD response of 
UF-JNJ. By projecting the continuum Hamiltonian onto 
a suitable set of localized states we reduced the problem 
to the solution of a discrete system in which the electro- 
magnetic field is described in terms of Peicrls phases. The 
latter provide the basic quantities to construct a density 
functional theory of superconducting (and as a special 
case normal) systems. We proved that under reasonable 
conditions the TD bond current and pairing density of an 
interacting system driven out of equilibrium by Peierls 
phases j(t) can be reproduced in a system of noninter- 
acting KS electrons under the influence of Peierls phases 
7'(t) and pairing field A'(t) and that 7'(t) and A'(t) are 
unique. We considered the KS system initially in equilib- 
rium at given temperature and chemical potential when 
at time t = an external electromagnetic field is switched 
on. To calculate the response of the system at times t > 
we used a non-equilibrium formalism in which the normal 
and anomalous propagators are defined on an extended 
Keldysh contour that includes a purely imaginary (ther- 
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FIG. 9: Photon-assisted MAR in a chain of 21 atomic sites. 
The equilibrium parameters are the same as in Fig. [5] An ac 
bias Ur = t/r sin(u> r £) in lead R is superimposed to a dc bias 
Ul = 0.8 in lead L. The panels show the contour plots of the 
bond-current J n ,n+i(t) across the atomic bonds of region C 
for different values of U t = 0.0, 0.1, 0.3, 0.5 and uj t = 0.4. 
The results are obtained with a time-step S = 0.05, cut-off 
A = 4 and a number of scattering states N p ,l = N p ,r = 1232. 



mal) path going from to —ifj. We showed that the 
solution of the equations of motion for the NEGF are 
equivalent to first solve the static BdG equations and 
then the TD BdG equations. It is worth emphasizing 
that in TDSCDFT the BdG equations do not follow from 
the BCS approximation and that their solution yields the 
exact bond-current and pairing density of an interacting 
system provided that the exact KS Peierls phases and 
pairing field are used. 

For systems consisting of Af superconducting leads in 
contact with a finite region C and driven out of equilib- 
rium by a longitudinal electric field a numerical algorithm 
is proposed. The initial eigenstates are obtained from a 
recent generalized wave-guide approach properly adapted 
to the superconducting case^ The initial states are prop- 
agated in time using an embedded Crank-Nicholson al- 
gorithm which is norm-conserving, accurate up to sec- 
ond order in the time-step and that exactly incorpo- 
rates transparent boundary conditions. The propagation 
scheme reduces to the one of Refs. I37lf38l in the case of 
normal leads. 

The method described in this work allows for obtain- 
ing the TD current across an UF-JNJ and hence to fol- 
low the time evolution of several AR until the Joscphson 
regime sets in. As a first calculation of these kind we 
explored in detail the popular single-level QD model in 
the weak and intermediate coupling regime. We demon- 
strated that the transient time increases with decreasing 
bias and provided a quantitative picture of the MAR. The 



rich structure of the transient regime is due to the ABS 
which play a crucial role in the relaxation process. For dc 
pulses we showed that ABS can be exploited to generate 
zero bias ac currents of tunable frequency. Furthermore, 
irradiating the biased system with a microwave field of 
appropriate frequency the ABS give rise to a long-living 
transient resonant effect. The transient regime increases 
also with the length of the junction. We considered one- 
dimensional atomic chains coupled to superconducting 
leads under dc and ac biases. Here we showed that in 
conditions of perfect AR there exists a clear-cut transient 
scenario for MAR. For biases Ul = 2A/n the dominant 
scattering channel is the n-th order AR and the transient 
regime lasts for about nN/vc where N is the length of 
the chain and vc the electron velocity at the Fermi level. 
Similar considerations apply to photon assisted MAR. A 
more careful analysis of the transient regime is beyond 
the scope of the present paper. However such analysis is 
of utmost importance if the ultimate goal of supercon- 
ducting nanoelectronics is to use these devices for ultra- 
fast operations. 

The TD properties presented in this work have been 
obtained using rather simple, yet so far unexplored, mod- 
els. A more sophisticated description of the Hamilto- 
nian is, however, needed for a quantitative parameter- 
free comparison with experiments. Theoretical advances 
also involve the development of approximate function- 
als for the self-consistent calculation of the TD pairing 
potential and Peierls phases. Self-consistent calculations 
have so far been restricted to equilibrium S-D-S mod- 
els with a point-like attractive interaction treated in the 
BCS approximation* 7 ^— For biased systems, however, 
the pairing potential and Peierls phases must be treated 
on equal footing and a first step in this direction would 
be the BCS approximation for the pairing field and the 
Hartree-Fock approximation for the Peierls phases. More 
difficult is the study of UF-JNJ in the Coulomb blockade 
regime for which electron correlations beyond Hartree- 
Fock must be incorporated. 

Finally, the approach presented in this work is not lim- 
ited to two terminal systems. The coupling of the cen- 
tral region to a third normal lead, or gate, allows for 
controlling the Joscphson current by varying the gate 
voltagei 25 ' 76 ' 77 These systems can be potentially used for 
fast switches and transistors ) 78 ' 79 and a microscopic un- 
derstanding of their ultrafast properties is therefore nec- 
essary to optimize their functionalities. 



Appendix A: Calculation of the embedding matrices 

Without loss of generality we include few layers of each 
lead in the explicitly propagated region C. Then, the 
embedding matrix Q*-™-* is zero everywhere except in the 
block of dimension 2N" cll x 2N" cll which is connected to 
the a lead. Denoting with q( m ) such non- vanishing block 
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in Q (m) we have 



Jm) 



m+1 



(Al) 



o.o 



where the subscript (0, 0) denotes the first diagonal block 
(supercell with j = 0) of the matrix in the square brack- 
ets. We notice that from Eq. (|75)l the matrix H aa is the 
same as the matrix H QQ (0) in Eq. (|63[) but with renor- 



malized diagonal blocks h. Q = fa Q — /i<r Q . In order to 
compute the q^ m )'s we introduce the generating matrix 
function 



ci 0,2/) = L 



1 



(A2) 



o.o 



which can also be expressed in terms of continued matrix 
fractions 



xl a + iySh a +y 2 S 2 t a - 



1 



' xla + iySh a + y 2 6 2 (i a (x, y) 

r 



?l a + iy6h a + y 2 6 2 t a -±- tt 
ti = t a p(x,y)tl, 



(A3) 



where the last step is an implicit definition of pj^x,y). 

The q( m ) 's are obtained from the generating matrix func- 
tion as 



Jm) 



1 



x=y=l 



d_ + d_ in 

dx dy 

= t a pi m Hl (A4) 

Using the identity + ^(^^(i,!/) = 0, 

we derive the following recursive scheme 



- <5 2 V (qW + 2q( fc - 1 > + q( fc - 2 ))p(™- fc ) 

fc=0 

(A5) 

with p( m ) = q( m ) = for m < 0. The above relation can 

—a —a 

be used to calculate q^"' provided that all p^ are known 
for k < m. To obtain p(°) we can use Eq. (|A3[) with 
x = y = 1 in which the continued fraction is truncated 
after a number -/Vievei of levels. Convergence can be easily 
checked by increasing Mevei- 

Appendix B: Calculation of the boundary term 



From Eq. ([81]) we see that in order to propagate an 
cigenstate of H — [ia_ we need to know the boundary 
term defined in Eq. ([52")). The state <E>(°) can be either a 
scattering state or an ABS. As shown in Section llll Al thc 
projection onto lead a of a generic cigenstate with energy 
E can be written as a linear combination of states of the 
form 



$t(m=s,j,a) = Z%(s)e 



ikj 



(Bl) 



where the amplitudes Z? satisfies the eigenvalue equation 



- M <rJ Zt = EZ%. (B2) 

In the following we show how to compute the action of 
the operator H Cct (0)g™ Q (l QQ + g ( 
the Nambu vector in region C 



on We define 



a(m) 
C,k = 



2H CQ (0) 



(1. 



m+1 k J 



(B3) 



from which the boundary term can easily be extracted by 
taking the appropriate linear combination of the ^'qT^ 
and then multiplying by —i5za n \ see Eq. (|82j) . Since re- 
gion C includes few layers of the leads the vector ™ is 
zero everywhere except for the components correspond- 
ing to orbitals in contact with lead a. If we call 
the vector with such components from Eq. (|B3[) we can 
write 



,a(m) 
IJ C,k 



2t„ 



m+1 



2« 



(B4) 

where the subscript j = in the square brackets denotes 
the vector of dimension 2N" cll with components given by 
the projection of the full vector onto the first (j = 0) 
supercell. As for the embedding matrices we introduce 
the generating function 



V k a (x,y) = 



1 



x\ r 



iySU c 



-n 



(B5) 



3=0 
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from which the V k 
tives 



v k — I 

K ml 



a(m) 



are obtained via multiple deriva- 



d_ d_ 

dx dy 



V k a (x,y) 



(B6) 



x=y=l 



The generating function can be obtained as follows. Tak- 
ing $^ as in Eq. (|B1[) and exploiting the property in Eq. 
(|B2|) it is easy to realize that 



H $f 



= (E- S^e-^tl) . , (B7) 



where the subscript j denotes the vector of dimension 
with components given by the projection of the 
full vector onto the j-th supercell. Then, multiplying the 
Dyson identity 



1 



1 

x 



iyS 



1 



on the right by using Eq. (|B7|) and solving for 
V^(x,y) we obtain the following result 

1 + iySe~ lk p (x, y)tj_ 
x + lyohj 



where p (x,y) is the generating function defined in Eq. 

(|A3|) . The quantity V^ m ^ can now be obtained from Eq. 
(IB6D and reads 



<m) (1 - iSEY 



_ 4fe ^ (i-iSEy 



v 7 n— 



(1 + iSE)™-^ 1 



x (Ei ra) +Pi" _1) )^ Q ^ B1 °) 
H QQ (B8) This conclude the calculation of the boundary term. 
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